Raport z analizy zmiany w rozmiarze wyławianego Śledzia Oceanicznego

Abstrakt

Wprowadzenie

Na przestrzeni ostatnich lat zauważono stopniowy spadek rozmiaru śledzia oceanicznego wyławianego w Europie. Do analizy zebrano pomiary śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Inicjalizacja środowiska

library(ggplot2)
library(dplyr)
library(DT)
library(GGally)
library(magick)
library(cowplot)

# ustawmy seed'a dla randoma w celu powtarzalności wyników
set.seed(1234)

Zainicjujmy również funkcje pomocnicze do prezentacji danych

prettyTable <- function(table_df, round_columns=numeric(), round_digits=2) {
    DT::datatable(table_df, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons",
                  options = list(dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel', 'pdf', 'print'))) %>%
    formatRound(round_columns, round_digits)
}

Import danych

Dane znajdują się w pliku sledzie.csv załączonym w repozytorium. Opisy ich atrybutów można znaleźć w treści zadania; dla wygody jednak zamieszczono go również poniżej:

Kolejne kolumny w zbiorze danych to:

• length: długość złowionego śledzia [cm];
• cfin1: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1];
• cfin2: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2];
• chel1: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1];
• chel2: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2];
• lcop1: dostępność planktonu [zagęszczenie widłonogów gat. 1];
• lcop2: dostępność planktonu [zagęszczenie widłonogów gat. 2];
• fbar: natężenie połowów w regionie [ułamek pozostawionego narybku];
• recr: roczny narybek [liczba śledzi];
• cumf: łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku];
• totaln: łączna liczba ryb złowionych w ramach połowu [liczba śledzi];
• sst: temperatura przy powierzchni wody [°C];
• sal: poziom zasolenia wody [Knudsen ppt];
• xmonth: miesiąc połowu [numer miesiąca];
• nao: oscylacja północnoatlantycka [mb].

Wiersze w zbiorze są uporządkowane chronologicznie.

sourceDataFileName <- "sledzie.csv"

Dokonajmy wstępnego wczytania danych - pierwszych 100 wierszy w celu automatycznego określenia klasy oraz podglądu danych

initialHerringData <- read.csv(sourceDataFileName, nrows = 100)
classes <- sapply(initialHerringData, class)
classes

        X    length     cfin1     cfin2     chel1     chel2     lcop1     lcop2      fbar      recr      cumf    totaln       sst       sal 
"integer" "numeric"  "factor"  "factor"  "factor"  "factor"  "factor"  "factor" "numeric" "integer" "numeric" "numeric"  "factor" "numeric" 
   xmonth       nao 
"integer" "numeric" 

Jak widać, większość danych została zakwalifikowana jako factor mimo, że z opisu wynika, że są to liczby. Wczytajmy zatem cały plik w odpowiednich klasach - wartości brakujące reprezentuje ?, natomiast jako separator dziesiętny przyjmujemy kropkę .. Plik nie zawiera również komentarzy. Szybka analiza pliku pozwoliła również na identyfikację, że zbiór składa się z około 52.6 tys wierszy - informację tą wykorzystamy w celu przyśpieszenia importu.

classes <- c("NULL", rep(c("numeric"), 13), "factor", "numeric")
herringData <- read.csv(sourceDataFileName, comment.char = "", header=TRUE, colClasses = classes, nrows=52600, dec = ".", na.strings = c("?"))

Omówienie danych

Sprawdzmy czy rozmiary zbioru pokrywają się z oczekiwaniami

dim(herringData)

[1] 52582    15

Wygląda na to, że dane zostały poprawnie wczytane. Zobaczmy je więc

prettyTable(head(herringData, n = 100))

Następnie przeprowadzmy krótkie podsumowanie

summary(herringData)

     length         cfin1             cfin2             chel1            chel2            lcop1              lcop2             fbar       
 Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000   Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
 1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469   1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
 Median :25.5   Median : 0.1111   Median : 0.7012   Median : 5.750   Median :21.673   Median :  7.0000   Median :24.859   Median :0.3320  
 Mean   :25.3   Mean   : 0.4458   Mean   : 2.0248   Mean   :10.006   Mean   :21.221   Mean   : 12.8108   Mean   :28.419   Mean   :0.3304  
 3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936   3rd Qu.:11.500   3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
 Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000   Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
                NA's   :1581      NA's   :1536      NA's   :1555     NA's   :1556     NA's   :1653       NA's   :1591                     
      recr              cumf             totaln             sst             sal            xmonth           nao          
 Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77   Min.   :35.40   8      : 9920   Min.   :-4.89000  
 1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.60   1st Qu.:35.51   10     : 7972   1st Qu.:-1.89000  
 Median : 421391   Median :0.23191   Median : 539558   Median :13.86   Median :35.51   7      : 6922   Median : 0.20000  
 Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87   Mean   :35.51   9      : 5714   Mean   :-0.09236  
 3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16   3rd Qu.:35.52   6      : 4218   3rd Qu.: 1.63000  
 Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73   Max.   :35.61   5      : 3736   Max.   : 5.08000  
                                                       NA's   :1584                    (Other):14100                     

Jak widać, śledzie mają długość w przedziale 19-32.5 cm, z medianą ustaloną dla wartości 25.5 cm, najwięcej śledzi wyławia się w sierpniu, październiku oraz lipcu. Można również stwierdzić, że liczba brakujących danych dla pól cfin, chel1, lcop oraz sst jest porównywalna; dla innych column wszystkie dane są określone. Możemy również zwrócić uwagę na wielkość połowów - ich liczba mieści się w przedziale 14.4 - 101.6 tys z medianą na poziomie 54 tys oraz średnią 51.5 tys. Rozkład ten jest więc prawdopodobnie zbliżony do normalnego, co z resztą możemy sprawdzić.

  ggplot(data=herringData, aes(x=totaln)) +
  geom_density() +
  labs(x = "Wielkość połowu", y = "Gęstość") +
  theme_bw()

Niekoniecznie - nie przypomina on żadnego z popularnych rozkładów.

Mamy za zadanie zweryfikować, czy przez ostatnie kilka lat rozmiar śledzia spadł. Przydatny może okazać sie więc wykres prezentujący średni rozmiar na rok. Jest jednak jeden problem - w danych mamy jedynie miesiąc połowu, informację o tym, że są to dane chronologicznie posortowane oraz pochodzą z ostatnich 60 lat. Należy więc w miarę możliwości dodać rok do danych bazując na tych informacjach. Sprawdzmy jednak czy nie ma wśród tych danych brakujących wartości

sum(is.na(herringData$xmonth))

[1] 0

brak - możemy więc przejść dalej. Dane są chronologiczne, także powinniśmy być w stanie wywnioskować rok na podstawie zmiany miesiąca na mniejszy

evaluateYearByMonthChange <- function(fishing) {
  curYear <- 0
  recentMonth <- 0
  l <- lapply(fishing, function(v) {
    v<- as.integer(v)
    if (v < recentMonth)
      curYear <<- curYear + 1
    recentMonth <<- v
    curYear
  })
  unlist(l)
}

herringDataWithYearByMonthChange <- mutate(herringData, year=evaluateYearByMonthChange(xmonth))
herringDataWithYearByMonthChange

Pogrupujmy je po roku dzięki funkcji

statsInYear <- function(dataWithYear) {
  dataWithYear %>%
    group_by(year) %>%
    summarise_at(vars(length), funs(mean(., na.rm=TRUE), min = mean - sd(., na.rm= TRUE), max = mean + sd(., na.rm= TRUE)))
}

meanLengthInYearChart <- function(dataWithYear) {
  meanPerYear <- statsInYear(dataWithYear)
  ggplot(data=meanPerYear, aes(x = year, y = mean)) +
    geom_ribbon(aes(ymin = min, ymax = max), alpha = 0.5, fill = "darkseagreen3", color = "transparent") +
      geom_line(color = "aquamarine4", lwd = 0.7) + 
    labs(x = "Rok", y = "Długość śledzia") +
  theme_bw()
}

meanLengthInYearChart(herringDataWithYearByMonthChange)

Założenie o chronologii jak widać okazało się błędne, przynajmniej w zakresie kolejności miesięcy. Jest również pole recr mówiące o połowie w roku. Widzimy, ze występuje pewna powtarzalność, więc jeśli dane byłyby chronologiczne, a zarazem ich rozpiętość jest całkiem spora, powinno to pozwolić na zgrupowanie.

evaluateYear <- function(fishing) {
  curYear <- 0
  recentYearFishing <- 0
  l <- lapply(fishing, function(v) {
    if (v != recentYearFishing)
      curYear <<- curYear + 1
    recentYearFishing <<- v
    curYear
  })
  unlist(l)
}

herringDataWithYearByRecr <- mutate(herringData, year=evaluateYear(recr))
herringDataWithYearByRecr

meanLengthInYearChart(herringDataWithYearByRecr)

Niestety, ponownie porażka. W takim razie może po prostu sprawdzmy czy występuje trend w czystych danych:

ggplot(herringData, aes(x=as.numeric(row.names(herringData)), y=herringData$length)) +
  geom_line(color = "aquamarine4") +
  labs(x = "Rok", y = "Długość śledzia") + 
  theme_bw()

Dobrze, jakąś zależność widać - może aby być jej pewnym uśrednijmy dane stosując prawo wielkich liczb oraz wiedzę o tym, że mamy do czynienia z danymi z okresu 60 lat

expectedYears <- 60
expectedRowsPerYears <- nrow(herringData) %/% 60
herringDataWithYearByPeriod <- herringData %>%
  mutate(year = 1:n() %/% expectedRowsPerYears)
herringDataWithYearByPeriod

meanLengthInYearChart(herringDataWithYearByPeriod)

Error in mean.default(length, na.rm = TRUE, function (x, na.rm = FALSE,  : 
  'trim' must be numeric of length one

Widzimy więc, że w rzeczywistości istnieje istotny spadek średniej długości śledzia na przełomie lat. Możemy to zwizualizować.

img <- image_read("herring.svg")
img <- image_trim(img)
heightRatio <- 110/475
maxWidth = 542
maxHeight <- heightRatio * maxWidth
herringsImgs <- statsInYear(herringDataWithYearByPeriod) %>%
  apply(1, function(r) {
    newWidth <- round((r[[2]] - 10) * 30)
    widthDif <- round((maxWidth - newWidth)/2)
    heightDif <- round(widthDif * heightRatio)
    image_scale(img, newWidth) %>%
        image_border("white", paste(c(widthDif, "x", heightDif), collapse = "")) %>%
        image_annotate(paste(c("Rok: ", r[[1]]), collapse = ""), color = "black", size = 30)
  })
frames <- image_morph(image_join(herringsImgs), frames = 1)
anim <- image_animate(frames, fps = 10)
image_write(anim, "herring.gif")

Zidentyfikujmy możliwe powody takiego zachowania wyznaczając macierz korelacji.

Wcześniej jednak należy zmienić typ kolumny xmonth z powrotem na liczbę i przygotować dane.

herringDataWithNumericMonth = mutate(herringData, xmonth = as.numeric(xmonth))

Zobaczmy wynik w bardziej przyjaznej formie

ggcorr(
  herringDataWithNumericMonth,
  nbreaks = 9,
  label = TRUE,
  label_size = 3,
  color = "grey50",
  layout.exp = 1,
  hjust = 0.75,
)

Widać dużą korelację pomiędzy parami chel1 - lcop1, chel2 - lcop2, fbar - cumf, cfin2 - lcop2 oraz cumf - totaln. Na samą długość w największym stopniu wpływa sst, nao oraz fbar. Widzimy że lcop1, lcop2 oraz cumf mogą zostać pominięte w procesie dalszej analizy jako redundantne. Teoretycznie moglibyśmy również usunąć miesiąc ze względu na powtarzalność, ale może mieć on wpływ na cykl rozwojowy, więc pozostawmy go póki co.

filteredHerringDataWithNumericMonth <- select(herringDataWithNumericMonth, -c(lcop1, lcop2, cumf))

Zobaczmy jeszcze jak przebiega zależność pomiędzy najbardziej obiecującymi zmiennymi na długość; kolorami wyróżniono również miesiące

numericChart <- function(colName) {
  ggplot(herringData, aes(x = length, y = herringData[[colName]], color=xmonth)) +
    geom_jitter(size = 1.5, stat = "identity") +
    labs(x = "Długość", y = colName)
}

plot_grid(
  numericChart("chel1"),
  numericChart("fbar"),
  numericChart("sst"),
  numericChart("nao")
)

Powiązanie rosnącego sst oraz nao jasno wpływa negatywnie na długość, przy fbar również można dostrzec pozywytną korelację.

Sprawdzmy również zależność pomiędzy miesiącem połowu a długością; w macierzy korelacji nie mogliśmy tego do konca zweryfikować, jednak ze względu na cykliczność było to również utrudnione (miesiąc w tym wypadku jest sztucznie liczbą).

numericChart("xmonth")

Wydaje się, że nie ma żadnej znaczącej zależności pomiędzy miesiącem połowu a długością śledzia